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Abstract 

An extension of the restricted Delaunay-refinement algorithm for surface mesh generation is described, where a new 
point-placement scheme is introduced to improve element quality in the presence of mesh size constraints. Specifically, 
it is shown that the use of off-centre Steiner points, positioned on the faces of the associated Voronoi diagram, typically 
leads to significant improvements in the shape- and size-quality of the resulting surface tessellations. The new algorithm 
can be viewed as a Frontal-Delaunay approach - a hybridisation of conventional Delaunay-refinement and advancing- 
front techniques in which new vertices are positioned to satisfy both element size and shape constraints. The performance 
of the new scheme is investigated experimentally via a series of comparative studies that contrast its performance with 
that of a typical Delaunay-refinement technique. It is shown that the new method inherits many of the best features 
of classical Delaunay-refinement and advancing-front type methods, leading to the construction of smooth, high quality 
surface triangulations with bounded radius-edge ratios and convergence guarantees. Experiments are conducted using a 
range of complex benchmarks, verifying the robustness and practical performance of the proposed scheme. 

Keywords: Surface mesh generation, Delaunay-refinement, Advancing-front. Frontal-Delaunay, Off-centre Steiner 
points. Element quality 


1. Introduction 

Surface mesh generation is a key component in a va¬ 
riety of computational modelling and simulation tasks, 
including many forms of computational engineering and 
numerical modelling, a range of problems in computer 
graphics and animation, and various data visualisation 
applications. Given a geometric domain described by a 
bounding surface E embedded in the surface trian¬ 
gulation problem consists of tessellating E into a mesh 
of non-overlapping triangular elements, such that all geo¬ 
metrical, topological and user-defined constraints are sat¬ 
isfied. While each use-case contributes its own set of spe¬ 
cific considerations, it is typical to require that surface 
tessellations: (i) consist of elements of high shape-quality, 
(ii) provide good geometrical and topological approxima¬ 
tions to the underlying surface E, and (iii) satisfy a set 
of user-specified element sizing constraints. While various 
strategies have been developed to solve the surface mesh¬ 
ing problem in the past, a new algorithm is developed in 
this study with the aim of improving the quality of the 
resulting triangulations. 
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1.1. Related Work 

Mesh generation is a broad and evolving area of re¬ 
search. A number of algorithms have been developed to 
tackle various meshing tasks, including grid-overlay meth¬ 
ods 013 ] in which a mesh is generated by warping and re¬ 
fining a semi-structured background grid, advancing-front 
techniques imaieiiT] which grow a mesh incrementally in 
a layer-wise fashion, and Delaunay-based strategies um 
uni [III na ns mmi m HD which incrementally refine an 
initially coarse Delaunay triangulation. 

The Delaunay triangulation [18] is imbued with a num¬ 
ber of attractive theoretical properties. In the context of 
mesh generation, for triangulations in M? and surface tes¬ 
sellations embedded in , it is known that such structures 
are topologically optimal m, maximising the worst-case 
element quality in the resulting mesh. The reader is re¬ 
ferred to [19] for additional details. As the name suggests, 
Delaunay-refinement schemes are based on the incremental 
refinement of an initially coarse bounding Delaunay tes¬ 
sellation. At each step of the algorithm, elements that vi¬ 
olate a set of geometrical, topological or user-defined con¬ 
straints are identified and the worst offending element is 
eliminated. This is achieved through the insertion of an 
additional Steiner vertex at the refinement point of the 
element - either the circumcentre of the element itself, or 
a point on an adjacent segment of the bounding geome¬ 
try. The original planar Delaunay-refinement methods of 
Chew [8] and Ruppert laiini have since been extended to 
support volumetric tessellations m and surface triangu¬ 
lations dSlITlITlI. 
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Delaunay-based methods have previously been applied 
to the surface meshing problem, via the so-called restricted 
Delaunay paradigm originally introduced by Edelsbrun- 
ner and Shah m- Given a surface S embedded in 
and a set of points X G S, the restricted Delaunay sur¬ 
face triangulation Del|s(X) is a sub-complex of the full¬ 
dimensional Delaunay tessellation Del(X), containing the 
set of 2-faces / G Del(X) that approximate the underly¬ 
ing surface. Boissonnat and Oudot ns HI] have shown 
that when the point-wise sampling X is sufficiently dense 
and well-distributed^ the restricted Delaunay triangulation 
Del |e(-^) is an accurate piecewise approximation of E, in¬ 
corporating several theoretical guarantees of fidelity. Given 
such behaviour, a Delaunay-refinement algorithm for sur¬ 
face mesh generation proceeds as per the planar case, with 
an initially coarse restricted Delaunay triangulation, in¬ 
duced by a sparse sampling X G S, being incrementally 
refined via the introduction of new Steiner vertices posi¬ 
tioned on the surface E. This process is discussed in detail 
in Section 121 

Advancing-front techniques are an alternative approach 
to mesh generation, in which elements are created one-by- 
one, starting from a set of frontal facets initialised on the 
bounding geometry. The mesh is marched inwards layer- 
by-layer, with new elements created by carefully position¬ 
ing vertices adjacent to facets in the frontal set. Following 
the placement of new elements, the set of frontal facets 
is updated to incorporate the new mesh topology. Ele¬ 
ments are added incrementally in this fashion until a com¬ 
plete tessellation of the domain is obtained. Such methods 
are applicable to surface meshing problems [51 0 [7], al¬ 
though several non-trivial issues of robustness are known 
to exist [7]. Fundamentally, advancing-front techniques 
are heuristic in nature, and do not typically offer theo¬ 
retical guarantees of convergence or bounds on element 
quality. Nonetheless, when such methods are successful, 
they are known to produce very high-quality tessellations. 

Due to their theoretical robustness, Delaunay-refine¬ 
ment schemes are a popular choice for practical meshing 
algorithms. On the other hand, it is well known that 
optimised advancing-front type methods can often pro¬ 
duce meshes with superior element quality characteris¬ 
tics, though their heuristics can sometimes break-down 
in practice. In this study, a new Frontal-Delaunay al¬ 
gorithm is presented which seeks to combine the bene¬ 
fits of classical Delaunay-refinement and advancing-front 
type approaches. While such strategies have been pursued 
previously in the case of planar and/or volumetric refine¬ 
ment algorithms, it is believed that the present study is 
the first application of such ideas to the restricted surface 
meshing problem directly. The new algorithm is designed 
to produce smooth, high-quality triangulations consistent 
with an advancing-front scheme, whilst also inheriting the 
theoretical robustness of Delaunay-based techniques. It 
is expected that this new algorithm may be of interest to 
users who place a high premium on mesh quality, including 
those operating in the areas of computational engineering 


and numerical simulation. 

A brief review of the traditional surface Delaunay-refine¬ 
ment scheme is presented in Section along with a dis¬ 
cussion of the underlying theoretical constructs, including 
the restricted Delaunay tessellation. The new Frontal- 
Delaunay algorithm is presented in Section focusing 
in-detail on the off-centre point-placement scheme pro¬ 
posed in this work. In Section a detailed comparison 
between the conventional Delaunay-refinement and pro¬ 
posed Frontal-Delaunay schemes is presented, contrasting 
output quality and computational performance. 

2. Restricted Delaunay-refinement Techniques 

Delaunay-refinement algorithms operate by incremen¬ 
tally adding new Steiner vertices to an initially coarse re¬ 
stricted Delaunay triangulation of the underlying surface. 
Contrary to typical planar algorithms [susiin], refinement 
schemes for surface meshing are designed not only to en¬ 
sure that the resulting mesh satisfies element shape and 
size constraints, but that the geometry and topology of 
the mesh itself is an accurate piecewise approximation of 
the underlying surface. 

Before discussing the details of Delaunay-refinement al¬ 
gorithms for surface meshing problems, a number of im¬ 
portant theoretical concepts are introduced: 

Definition 1. Let F be a smooth surface embedded in 
enclosing a volume C Let Del(X) be a full¬ 
dimensional Delaunay tessellation of a point-wise sample 
X C F and Vor(X) be the Voronoi diagram associated 
with Del(X). The restricted Delaunay surface triangula¬ 
tion Del|E(X) is a sub-complex of Del(X) including any 
2-face / G Del(X) associated with an edge vj G Vor(X) 
such that vj DF 7 ^ 0. The restricted Delaunay volume tri¬ 
angulation Del |o(X) is a sub-complex of Del(X) including 
any 3-simplex r G Del(X) associated with an internal cir- 
cumcentre c G D. 

It has been shown by a number of authors, including 
Amenta and Bern [ 21 ], Cheng, Dey, Levine [ 22 ], Cheng, 
Dey and Ramos na, Cheng, Dey, Ramos and Ray [23], 
Boissonnat and Oudot nsms] and Cheng, Dey and Shewchuk 
[19], that given a sufficiently dense point-wise sampling 
of the surface, X G F, the restricted Delaunay tessella¬ 
tion Del |e(-^) is guaranteed to be both geometrically and 
topologically representative of the underlying surface F. 
Specifically, it is known that, when the restricted tessel¬ 
lation is a so-called loose e-sample [El US], Del|E(X) is 
homeomorphic to the surface F, the Hausdorff distance 
H(Del |e(X), F) is small and that Del|E(X) provides a 
good piecewise approximation of the geometrical prop¬ 
erties of F, including its normals, area and curvature. 
Similar theoretical guarantees extend to the associated re¬ 
stricted volume tessellation Del|o(X). The properties of 
restricted Delaunay tessellations are well documented in 
the literature, and a full account is not given here. The 
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Figure 1: Restricted tessellations for a curved domain in showing (i) the curve S and enclosed area (ii) the Delaunay tessellation 
Del(X) and Voronoi diagram Vor(X), and (iii) the restricted curve and area tessellations Del| 2 (-^) and Del ln(^)- 





reader is referred to the detailed expositions presented in 
[Tang and m for further details and proofs. 

Definition 2. Let Del |e(-^) be a restricted Delaunay tri¬ 
angulation of a smooth surface S. Given a 2-simplex 
/ G Del IE 5 any circumball of / centred on the sur¬ 
face H is known as a surface Delaunay hall of /, denoted 
SDB(/). 

Surface Delaunay balls are centred at the intersection 
of the Voronoi diagram Vor(X) with the surface S. Specif¬ 
ically, each 2-face / G Del |e(-^) is associated with an or¬ 
thogonal edge in the Voronoi diagram vj G Vor(V), which 
is guaranteed, by definition, to intersect with the surface S 
at least once. In some cases, especially when the sampling 
V G S is relatively coarse, there may be multiple surface 
intersections associated with a given 2-face / G Del |e(-^)- 
In such cases, it is typical to consider only the largest as¬ 
sociated surface ball. See Figure for an illustration of 
the surface Delaunay ball for a general facet /. 

Definition 3. Let Del |e(-^) be a restricted Delaunay sur¬ 
face triangulation of a smooth surface S. Given a 2- 
simplex / G Del|E(V), the surface discretisation error 
e(/) is the Euclidean distance between the centres of the 
largest surface Delaunay ball of / and its diametric ball. 

The surface discretisation error is a measure of how 
well the restricted triangulation Del | e approximates 

the underlying surface E geometrically. Gonsidering a 2- 
face / G Del |e(-^), the surface discretisation error e(/) is 
a measure of the distance between the face / and the fur¬ 
thest adjacent point on the surface E. This measure can 
be thought of as a one-sided discrete Hausdorff distance, 
e(/) = H(Del |e(V), E), defined from a set of representa¬ 
tive points on the triangulation Del |e(V) to the surface E. 
Glearly, if the triangulation Del |e(-^) is a good piecewise 
representation of the surface E, the surface discretisation 
error e(/) should be small. See Figure]^ for a description 
of the surface discretisation error for a facet / G Del |e(-^)- 

Definition 4. Given a d-simplex r, its radius-edge ratio, 
p(r), is given by: 

P(t) = -R/l|eo||, (1) 


where R is the radius of the diametric ball of r and ||eo|| 
is the length of its shortest edge. 

The radius-edge ratio is a measure of element shape- 
quality. It achieves a minimum, p(r) = i/x/s, for equilat¬ 
eral triangles and increases toward Too as elements tend 
toward degeneracy. For 2-simplexes, the radius-edge ratio 
is robust and can be related to the minimum plane angle 
^min between adjacent edges, such that p(t) = ^(sin(^ min)) ^ 
Due to the summation of angles in a triangle, given a min¬ 
imum angle Omin the largest angle 6 >max is clearly bounded, 
such that 6 >max < 180° — 26>min- 

2.1. An Existing Delaunay-refinement Algorithm 

The development of provahly-good Delaunay-refinement 
schemes for surface-based mesh generation is an ongoing 
area of research. An algorithm for the meshing of closed 
2 -manifolds embedded in adapted largely from the 
methods presented by Boissonnat and Oudot in [nmni 
is presented here. This method is largely equivalent to 
the CGALMESH algorithm, available as part of the CGAL 
package, and summarised by Jamin, Alliez, Yvinec and 
Boissonnat in m- A similar algorithm is also described 
by Gheng, Dey and Shewchuk in m and Dey and Levine 
in [24], while a variant developed for the reconstruction 
of medical images is presented by Foteinos, Ghernikov and 
Ghrisochoides in [25] . The algorithm presented in this sec¬ 
tion is referred to as the ‘conventional’ Delaunay-refinement 
approach throughout, due to its direct use of circumcentre- 
based Steiner vertices. 

As per Jamin et al. m, the development of the conven¬ 
tional Delaunay-refinement algorithm is geometry-agnostic^ 
being independent of the specific description used for the 
underlying surface E. It is required only that the geome¬ 
try support a so-called oracle predicate that can be used to 
compute the intersection of a given line segment with the 
surface E. The Frontal-Delaunay algorithm presented in 
Section]^ additionally requires that the oracle compute in¬ 
tersections between a disk of a given radius and the surface 
E. Such constructions will be discussed in further detail in 
Section]^ While a very general class of surface description 
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Figure 2: The surface Delaunay ball for a restricted 2-face / G Del showing (i) placement of the surface ball at the intersection of the 

dual edge vj G Vor(X) and the surface S, and (ii) associated radius r(/) and surface discretisation error e(/). 



is supported at the theoretical level, in this study, atten¬ 
tion is restricted to the development of so-called remeshing 
operations, in which the underlying surfaces H are speci¬ 
fied as existing 2-manifold triangular complexes V. This 
restriction is made in the present study for convenience 
only, and future work is intended to focus on more gen¬ 
eral surface descriptions, including implicit and analytic 
functions, in addition to those that contain sharp creases. 

Following Jamin et al. m, the Delaunay-refinement 
algorithm takes as input a surface domain, described by 
a smooth 2-manifold S, an upper bound on the allowable 
element radius-edge ratio p, a mesh size function h(x) de¬ 
fined at all points on the surface S and an upper bound on 
the allowable surface discretisation error e(x). The input 
surface H encloses a bounded volume ft. The algorithm 
returns a triangulation T|e of the surface S, where T|e 
is a restricted Delaunay surface triangulation of a point- 
wise sampling X G S, such that T|e = Del|E(X). As a 
by-product, the algorithm also returns a coarse triangu¬ 
lation T\n of the enclosed volume ft, where T\n is a re¬ 
stricted Delaunay volume triangulation T\n = Del |o(X). 
Both Del|^(X) and Del|f^(X) are sub-complexes of the 
full-dimensional Delaunay tessellation Del(X). Note that 
Del|E(X) is a triangular complex, while Del|f^(X) and 
Del(X) are tetrahedral complexes. The method is sum¬ 
marised in Algorithm |2.1 1 

The Delaunay-refinement algorithm guarantees, firstly, 
that all elements in the output surface triangulation / G 
T|e satisfy element shape constraints, p{f) < p, element 
size constraints h{f) < h{xf) and surface discretisation 
bounds e(/) < e(x/). Furthermore, for sufficiently small 
mesh size functions h(x), the surface triangulation T|e 
is guaranteed to be a good piecewise approximation of 
the underlying surface S, exhibiting both geometrical and 
topological convergence as h(x) ^ 0, consistent with the 
characteristics of restricted Delaunay tessellations outlined 
in Definition [T] 

The Delaunay-refinement algorithm begins by creating 
an initial point-wise sampling of the surface Xq G S. Ex¬ 
ploiting the discrete representation available for E, in this 
study the initial sampling is obtained as a well-distributed 
subset of the existing vertices Y eV, where V is the poly¬ 



hedral representation of the surface E. Care is taken to 
ensure that each connected component in V is sampled. 
Alternative initialisation methods making use of robust 
persistent-triangles techniques are described in [ms]. 

In the next step, the initial triangulation objects are 
formed. In the present work, the full-dimensional De¬ 
launay tessellation, Del(X), is built using an incremental 
Delaunay triangulation algorithm, based on the Bowyer- 
Watson technique [26]. The restricted surface and vol¬ 
umetric triangulations, Del|s(X) and Del|o(X), are de¬ 
rived from Del(X) by explicitly testing for intersections 
between edges in the associated Voronoi diagram Vor(X) 
and the surface E. These queries are computed efficiently 
by storing the surface definition V in a spatial-tree. 

The main loop of the algorithm proceeds to incremen¬ 
tally refine any 2-faces / G Del|^(X) that violate either 
the radius-edge, element size or surface discretisation re¬ 
quirements. The refinement process is priority scheduled, 
with triangles / G Del|E(X) ordered according to their 
radius-edge ratios p(/), ensuring that the element with 
the worst ratio is refined at each iteration. Individual el¬ 
ements are refined based on their surface Delaunay balls, 
with a triangle / G Del(X) eliminated by inserting the 
centre of the largest ball B(c, r)^^^ = SDB(/) into the tes¬ 
sellation Del(X). This process is a direct generalisation of 
the circumcentre-based insertion method associated with 
Ruppert’s algorithm for planar domains laiio]. 

As a consequence of changes to the full-dimensional 
tessellation following the insertion of a new Steiner ver¬ 
tex, corresponding updates to the restricted triangulations 
Del|E(X) and Del|o(X) are instigated, ensuring that all 
tessellation objects remain valid throughout the refine¬ 
ment process. The Delaunay-refinement algorithm termi¬ 
nates when all 2-faces / G Del|E(X) satisfy all radius- 
edge, size and surface discretisation thresholds, such thalQ 
p{f) < p, h{f) < ah{xf) and e(/) < e(x/), respectively, 
where the element size h{f) is proportional to the radius 
of the associated surface ball B(c,r)^^^ = SDB(/), such 


^The coefficient a = 4/3, ensuring that the mean element size does 
not, on average, undershoot the desired target size h(xj^). 
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thal[^/i(/) = v^r, and the target size h{'Kf) is sampled at 
the centre of the surface ball xj = c. 

3. Restricted Prontal-Delaunay Methods 

Frontal-Delaunay algorithms are a hybridisation of adv¬ 
ancing-front and Delaunay-refinement techniques, in which 
a Delaunay triangulation is used to define the topology 
of a mesh while new Steiner vertices are inserted in a 
manner consistent with advancing-front methodologies. In 
practice, such techniques have been observed to produce 
very high-quality meshes, inheriting the smooth, semi- 
structured vertex placement of pure advancing-front meth¬ 
ods and the optimal mesh topology of Delaunay-based ap¬ 
proaches. While Frontal-Delaunay methods have previ¬ 
ously been used by a range of authors in the context of pla¬ 
nar, volumetric and parametric surface meshing, includ¬ 
ing, for example studies by Ungor and Erten Rebay 
[28] , Mavriplis [29] , Frey, Borouchaki and George [30] , and 
Remade, Henrotte, Carrier-Baudouin, Bechet, Marchan- 
dise, Geuzaine and Mont on m, the authors are not aware 
of any previous investigations describing the application of 
such techniques to the surface meshing problem directly. 
The conventional advancing-front method, on the other 
hand, has been generalised to support surface meshing, 
as, for example, outlined in studies by Rypl [sa E] and 
Schreiner, Scheidegger, Fleishman and Silva [71 [33]. 

In these previous studies, the conventional planar adva¬ 
ncing-front methodology is extended directly to facilitate 
surface operations. Meshing proceeds via the incremen¬ 
tal introduction of a well-distributed set of vertices X 
positioned on the surface X G E. It should be noted 
that, contrary to the restricted Delaunay techniques in¬ 
troduced in previous sections, advancing-front methods 
typically maintain a partial, and possibly non-Delaunay 
manifold triangular complex T|e only - they do not con¬ 
struct a full-dimensional tessellation T|o. It should also 
be noted that, in addition to the usual limitations associ¬ 
ated with advancing-front strategies, serious issues of ro¬ 
bustness often afflict these techniques in practice due to 
the difflculties associated with the reliable evaluation of 
the requisite geometric intersection and overlap predicates 
for sets of non-planar elements. Additionally, for highly 
curved and/or poorly separated surface definitions it is 
known to be difficult to ensure the topological correctness 
of the output mesh T|e- Schreiner et al. [71|33] introduce 
a number of heuristic techniques in an effort to overcome 
these difflculties. 

3.1. Off-centres 

The new Frontal-Delaunay algorithm is based primar¬ 
ily on a generalisation of ideas introduced by Rebay, who. 


^The coefficient represents the mapping between the edge 
length and diametric ball radius for an equilateral element. Such 
scaling ensures that size constraints are applied with respect to mean 
edge length. 


in [28] , developed a planar Frontal-Delaunay algorithm in 
which new vertices are positioned along edges of the as¬ 
sociated Voronoi diagram. Rebay showed that new ver¬ 
tices can be positioned on Vor(X) according to an a priori 
mesh size function h(x), a strategy broadly consistent with 
conventional advancing-front techniques. While his algo¬ 
rithm maintains a Delaunay triangulation T = Del(X) 
of the current vertex set X, it is still fundamentally an 
advancing-front scheme - bereft of guarantees on the ele¬ 
ment shape quality p{r). Rebay reported that his scheme 
produced very high-quality output, typically outperform¬ 
ing conventional Delaunay-refinement in practice. 

Ungor and Erten have approached the problem from 
the flip-side, introducing the notion of generalised Steiner 
vertices for planar Delaunay-refinement methods. Their 
‘off-centre’ vertices are points lying along edges in the as¬ 
sociated Voronoi diagram, as per Rebay. In [34], Ungor 
showed that when refining a triangle r, its off-centre should 
be positioned, if possible, such that the new triangle a ad¬ 
jacent to the shortest edge in r satisfies the desired shape- 
quality target, such that p(cr) < p. Importantly, Ungor 
demonstrated that such a strategy typically leads to an 
improvement in the performance of Delaunay-refinement 
in practice, reducing the size of the output \T\. Ungor 
extended the guarantees derived for Ruppert’s Delaunay- 
refinement algorithm to his off-centre technique and showed 
that such bounds are typically improved upon in practice. 
Off-centre refinement is based on Ruppert’s Delaunay-refine¬ 
ment framework directly, involving modifications to the 
procedure used to refine low-quality triangles only. 

The use of generalised Delaunay-refinement strategies 
has also been explored by Ghernikov, Chrisochoides and 
Eoteinos in [35] [36] , in which Steiner points are positioned 
within a set of so-called selection-balls adjacent to ele¬ 
ment circumcentres. In [35], Ghernikov and Chrisochoides 
show that a family of ‘provably-good’ generalised two- 
and three-dimensional Delaunay-refinement schemes ex¬ 
ist, and can be realised via the specification of appropri¬ 
ate selection-ball radii parameters. It is unclear whether 
this selection-ball formalism incorporates the full set of 
Voronoi-type off-centres currently in use, including those 
introduced by Rebay [28], and Erten and Ungor m, in ad¬ 
dition to the strategies defined in the present study. Such 
a determination may form the basis for future work. 

3.2. Point-placement Preliminaries 

In this study, a generalisation of the ideas introduced 
by Rebay and Ungor is formulated for the surface mesh¬ 
ing problem - using off-centre Steiner vertices to simulate 
the vertex placement strategy of a conventional advancing- 
front approach, while also preserving the framework of 
a restricted Delaunay-refinement technique. The aim of 
such a strategy is to recover the high element qualities 
and smooth, semi-structured point-placement generated 
by frontal methods, while inheriting the theoretical guar¬ 
antees of Delaunay-refinement methods. Advancing-front 
algorithms typically incorporate a mesh size function h(x). 
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Algorithm 2.1 Surface mesh generation by restricted Delaunay-refinement 

1 

function DelaunaySurface(E, p, e(x), ^(x), T e, T n) 

1 

function RefineSimplex(/) > {surface refinement} 

2 

Form an initial sampling X G S such that X 

2 

Call SurfaceDelaunayBall(/, B(c, 


is well-distributed on E. 

3 

Form new Steiner vertex p about B(c,r)^^^. 

3 

Form Delaunay tessellation Del(X). 

4 

Insert Steiner vertex X ^ p, update Del(X) ^ X. 

4 

Form restricted objects Del s(-^) and Del n(^). 

5 

end function 

5 

Enqueue all restricted 2-simplexes Q|e ^ 

1 

function SurfaceDelaunayBall(/) > {surface ball} 


/ G Del a 2-simplex / is enqueued 

2 

Form Voronoi edge 'Vf orthogonal to 2-simplex /. 


if BadSimplex(/) returns true. 

3 

Form the set of associated surface Delaunay 

6 

while (Q s / 0) do > {main refinement sweeps} 


balls B(c, r)- for the restricted 2-simplex / G 

7 

Call RefineSimplex(/ ^ Q|e) 


Del e(A). Balls are centred about the set of 

8 

Update the restricted Delaunay tessel¬ 


surface intersections vj n E / 0. 


lations Del e(A) and Del q(X). 

4 

return B(c,r)^g^^, where rmax is maximal. 

9 

Update Q E to reflect changes to Del e(A). 

5 

end function 

10 

end while 

1 

function BadSimplex(/) > {termination criteria} 

11 

return TIe = Del|E(A) and r\n = Del|n(X) 

2 

return p{f) > p or e{f) > e(xp) or h{f) > h(x/) 

12 

end function 

3 

end function 


a function / : ^ M+ defined over the domain to be 

meshed, where h(x) represents the desired edge length ||e|| 
at any point x G S. This mesh size function typically in¬ 
corporates size constraints dictated by both the user and 
the geometry of the domain to be meshed. The construc¬ 
tion of appropriate mesh size functions will be discussed 
in subsequent sections, but for now, it suffices to note that 
h(x) is a ^f-Lipschitz function defined at all points on S 
with 0 < ^ < 1. 

The proposed Frontal-Delaunay algorithm is an ex¬ 
tension of the restricted Delaunay-refinement algorithm 
presented in Section modified to use off-centre rather 
than circumcentre-based refinement strategies. The basic 
framework is consistent with the Delaunay-refinement al¬ 
gorithm described previously, in which an initially coarse 
restricted Delaunay triangulation of a surface S is refined 
through the introduction of additional Steiner vertices X G 
S until all constraints are satisfied. A restricted surface 
triangulation T|e = Del|s(-^) is constructed as a sub¬ 
complex of the full-dimensional tessellation Del(X) and a 
coarse restricted volumetric tessellation T\n = Del|f^(X) 
is also available as a by-product. The constraints satis¬ 
fied by the Frontal-Delaunay algorithm are identical to 
those described previously for the restricted Delaunay- 
refinement scheme, with upper bounds on the radius-edge 
ratio p, surface discretisation error e(xj) and element size 
h(xj) all required to be satisfied for convergence. See Al¬ 
gorithm for a detailed outline of the method. 

3.3. Point-placement Strategy 

Given a surface facet / G Del|s(-^) marked for re¬ 
finement, the new Steiner vertex introduced to eliminate 
/ is an off-centre, constructed to adhere to local size and 
shape constraints. The off-centres introduced in this study 
involve the placement of two distinct kinds of Steiner ver¬ 
tices. Type I vertices, are equivalent to conventional 
element circumcentres, and are used to satisfy constraints 
on the element radius-edge ratios. Type II vertices, c^‘^\ 
are size-optimal points, designed to satisfy element siz¬ 
ing constraints in a locally optimal fashion. Adopting the 


generalised off-centre framework introduced by Ungor, the 
‘ideal’ location of the size-optimal off-centre is based 
on a consideration of the isosceles triangle cf formed about 
the short ‘frontal’ edge eo G /. The point is positioned 
to ensure that a satisfies local size constraints. 

Given a refinable 2-simplex / G Del|E(-^), the size- 
optimal Type II vertex is positioned at an intersection 
of the surface S and a plane V, where V is aligned with 
the local face of the Voronoi diagram Vor(X) associated 
with the frontal edge eo G /. The plane is positioned 
such that it passes through three local points on Vor(X): 
the midpoint of the frontal edge eo G /, the centre of the 
diametric ball of / and the centre of the surface Delaunay 
ball B(c,r)^^^ = SDB(/). The vertex is positioned 
such that it forms an isosceles triangle candidate a about 
the frontal edge eo, such that its size h{a) satisfies local 
constraints. Specifically, the altitude of a is computed 
from local mesh-size information, such that 

Ga = TXim{{hl - \\leof)i ,^/2ha) (2) 

^<T = 5(^(mi) + ft(m2)) (3) 

where the m^’s are the edge midpoints and is the 

altitude for an ‘ideal’ element. 

The position of the point is calculated by comput¬ 
ing the intersection of the surface S with a circle of radius 
Ucr, centred at the midpoint of the frontal edge eo G / 
and inscribed on the plane V. In the case of multiple in¬ 
tersections, the candidate point of closest alignment 
to the frontal direction vector d/ is selected. Specifi¬ 
cally, the point that maximises the scalar product 

— mo) • d/ is chosen, where mo is the midpoint of 
the short edge eo G / and the frontal direction vector df 
is taken from the midpoint mo to the centre of the surface 
ballB(c,r)^,, = SDB(/). 

For non-uniform h(x), expressions for the position of 
the point are weakly non-linear, with the altitude Oa- 
depending on an evaluation of the mesh size function at the 
edge midpoints h(m^) and visa-versa. In practice, since 
h(x) is guaranteed to be Lipschitz smooth, a simple it- 
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Figure 3: Placement of off-centre vertices, showing (i) the surface ball associated with a facet / G Del| 5 ](-^), (h) the plane V aligned with 
the local facet of Vor(X) associated with the frontal edge eo G /, (hi) placement of the size-optimal vertex 



erative predictor-corrector procedure is sufficient to solve 
these expressions approximately. The positioning of size- 
optimal Type II Steiner vertices is illustrated in Figure 

Using the size-optimal Type II point and the Type I 
point the final position of the refinement point c for 
the facet / is calculated. The point c is selected to satisfy 
the limiting local constraints, setting 

if and > ^||eo||), 

otherwise 

where — mo|| and — mo|| are 

distances from the midpoint of the frontal edge eo to the 
Type I and Type II vertices, respectively. 

The selection criteria is designed to ensure that the 
refinement scheme smoothly degenerates to that of a con¬ 
ventional circumcentre-based Delaunay-refinement strat¬ 
egy in limiting cases, while using a locally shape-optimal 
approach where possible. 

Specifically, the condition d*^^^ < d*^^^ guarantees that 
c lies no further from the frontal edge eo than the cen¬ 
tre of the surface Delaunay ball SDB(/). Such behaviour 
ensures that a conventional circumcentre-based scheme is 
selected when the element size becomes sufficiently small. 
Additionally, the condition d*^^^ > ^||eo|| ensures that the 
size-optimal scheme is only chosen when the edge eo is suf¬ 
ficiently small compared to the local mesh size function. 
This behaviour ensures that eo is a good ‘frontal’ edge 
candidate. 

3.4‘ Refinement Order 

In addition to the use of ‘off-centre’ point-placement 
schemes, the Frontal-Delaunay algorithm also introduces 
changes to the order in which elements are refined. To 
better mimic the behaviour of an advancing-front type 
method, elements are refined only if they are adjacent to 
an existing ‘frontal’ entity. In the case of surface facets / G 
Del the short frontal edge eo G / must be shared by 

at least one adjacent facet fj G Del|s(-^) that has ‘con¬ 
verged’ - satisfying its associated radius-edge, mesh-size 
and surface-error constraints. The idea of defining ‘frontal’ 
entities as a dynamic boundary between converged and 


un-converged elements in the mesh is a common feature 
of Frontal-Delaunay algorithms, with similar approaches 
used in [28l [29l [30l [31] . 

3.5. Theoretieal Guarantees 

Boissonnat and Oudot [151 US] have previously shown 
that the conventional restricted Delaunay-refinement al- 
gorithim presented in Section is guaranteed to: (i) ter¬ 
minate in a finite number of steps, (ii) produce elements 
of bounded radius-edge ratios p{f) < p, (iii) satisfy non- 
uniform user-defined mesh-sizing and surface error con¬ 
straints h{f) < h{'Kf)^ e(/) < e(xj), and (iv) provide 
a good geometrical and topological approximation to the 
underlying surface H when the mesh size function h(x) is 
sufficiently small. 

The behaviour of the new Frontal-Delaunay algorithm 
can be assessed using a similar framework. 

3.6. Termination & Convergenee 

Firstly, the termination and convergence of the algo¬ 
rithm is analysed: 

Proposition 1. Given a smooth surfaee S, a radius-edge 
threshold p > 1 and positive mesh-size and surfaee-error 
funetions h(x) > ho, e(x) > eo, ho,eo G the Frontal- 
Delaunay algorithm is guaranteed to terminate in a finite 
number of steps. 

Before seeking to prove Proposition itself, a number 
of intermediate results are first established. 

Lemma 1. Let D C be a bounded domain. Given a 
subset X C D, satisfying ||u—v|| > j for all pairs u,\- G X 
and some sealar 7 G the size of X is bounded, sueh 
that \X\ < jii for some eonstant p G Z+. 

This so-called paeking-lemma states that any point- 
set X in a bounded domain for which there exists a 
positive minimum separation distance between the points 
in X must, consequently, be finite. As such, if it can be 
shown that a meshing algorithm preserves such a minimum 
separation length between its vertices. Lemma can be 
invoked to prove that such a vertex-set is finite, and, as 
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a result, that termination of the algorithm is guaranteed. 
Lemma is stated here without proof - interested readers 
are referred to m for full details. 

Before examining the behaviour of the Frontal-Delaunay 
algorithm as a whole, the impact of a single refinement 
operation is analysed. Firstly, the behaviour of the vari¬ 
ous off-centre point-placement schemes introduced in Sec¬ 
tion [3]3] is examined: 

Proposition 2. Let the threshold on radius-edge ratios he 
p > 1 and Del|s(-^/c) he the restrieted surfaee triangula¬ 
tion indueed after k refinement operations. Given a low- 
quality surfaee faeet f G Del |e(-^/c); for whieh p{f) > p, 
the minimum edge length ||eo|| G Del|E(-^/c) 'is not de¬ 
creased following the insertion of a new Type I Steiner 
vertex at the eentre o/SDB(/). 

Proof Recalling that the surface facet / is locally Delau¬ 
nay, the insertion of a new Type I vertex Ck at the centre of 
the associated surface Delaunay ball B(c/c,r/c) = SDB(/) 
is guaranteed to create edges no shorter than the radius of 
the ball of /: 


||eo||/c+i > Tk- 


(5) 


Here ||eo||/c+i is the length of the shortest edge created 
by the insertion of the new vertex c/^. Dividing by the 
current minimum edge length, and using the definition of 
the radius-edge ratio: 


l|eo||fc+i Tk 

lleolk - lleolk 


(6) 


Clearly, when p > 1 the existing minimum length ||eo||/c is 
preserved by the insertion of the new vertex c/^. □ 

Proposition 3. Let h(x) > ho, ho G he a positive 
mesh-size funetion and Del|s(X/c) he the restrieted sur¬ 
faee triangulation indueed after k refinement operations. 
Given a refinahle surfaee faeet f G Del |e(-^/c) '^se of the 
Type II point-plaeement seheme is ‘deelined’ when f is suf- 
fieiently small. Speeifieally, seleetion of the Type II point- 
plaeement seheme requires that r^ > H, where r^ is the 
radius of the surfaee hall assoeiated with the faeet f and 
H is the length of the new frontal edge eandidates defined 
via expressions 

Proof. Given a surface facet /, the off-centre selection cri¬ 
teria expressed in @ requires that the Type II off-centre 
lies along a ‘safe’ region of the Voronoi diagram - 
bounded by the centre of the ball B{ck,rk) = SDB(/) 
such that: 


- moll < ||cfe - moll (7) 

where mo is the midpoint of the base edge ||eo||. Express¬ 
ing 0 in terms of the target edge lengths and surface ball 
radii: 


^i72_||leo||2<^,2_||l3^||2 


where H is the desired edge-length, computed from lo¬ 
cal mesh-size constraints 0 0 . This expressions can be 
further simplified, leading to a simple inequality for the 
selection of Type II vertices: 

H <rk ( 9 ) 

It can be seen that the Type II point-placement scheme 
is deelined when the surface ball associated with a given 
facet / is sufficiently small compared to the local target 
edge-length H. □ 

Using Propositions and the off-centre point-place¬ 
ment strategies utilised in the Frontal-Delaunay algorithm 
are shown to constitute a ‘safe’ set of refinement opera¬ 
tions, with the insertion of Type I vertices leading to a non¬ 
decreasing minimum edge-length and use of the Type II 
scheme declined once the mesh is sufficiently well refined. 
Using these results, termination of the full Frontal-Del¬ 
aunay algorithm can be re-visited: 

Proof of Proposition [TJ The Frontal-Delaunay algo¬ 
rithm refines any surface facet / G Del |e(-^) T- (i) h is of 
poor shape quality, such that p{f) > p, (ii) it is too large, 
such that h{f) > ah{'Kf), or (hi) it violates the local sur¬ 
face discretisation error threshold, such that e(/) > e(xj). 

(i) Using Proposition]^ it is known that the introduc¬ 
tion of Type II off-centres is declined once the mesh 
is sufficiently well refined. Type I points alone can 
therefore be relied upon to satisfy element shape- 
constraints. 

Given a radius-edge threshold p > 1, Proposition]^ 
states that refinement does not lead to a reduction 
in minimum edge length. Shape-based refinement 
therefore preserves the existing minimal edge length 
||eo ||/c, associated with either the initial sampling Xo G 
E, or the insertion of some previous vertex x/^, G X 
due to size- or surface-error-driven refinement. 

(ii) Noting that h(x) is positive, with h(x) > ho for some 
ho G M+, size-driven refinement is clearly deelined 
once the local element size is sufficiently small. Re¬ 
calling that h{f) = \/3r, where r is the radius of 
the surface ball associated with a given facet /, the 
mesh-sizing constraints h{f) < ah(xj) are satisfied 
once all r < («/U 3 )h(x/). 

(hi) Noting that e(x) is positive, with e(x) > eo for some 
eo G it is clear that surface-error driven refine¬ 
ment is deelined once the local element size is suffi¬ 
ciently small. Recalling that e(/) is the orthogonal 
distance between a given facet / and the centre of 
its associated surface ball B(c,r) = SDB(/), it is 
clear that e(/) < r, as SDB(/) must circumscribe 
the facet. In the worst-case, surface-error driven re¬ 
finement is declined once all r < e(x/). 


(8) 
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Using (i)-(iii), the progression of the Frontal-Delaunay 
algorithm can be examined. Starting with an initial ver¬ 
tex distribution Xq G S, both Type I and Type II point- 
placement schemes are used to satisfy mesh-size and surface- 
error constraints - with refinement continuing until all 
surface balls are sufficiently small, as per (ii)-(iii). Any 
remaining low-quality elements in violation of the radius- 
edge constraints are subsequently refined through the in¬ 
sertion of Type I vertices alone. This final phase is guar¬ 
anteed to preserve minimum edge-length. 

Using Lemma it is clear that the Frontal-Delaunay 
algorithm is guaranteed to produce a point-wise sampling 
X G S of finite size and must terminate in a finite number 
of steps as a consequence. 

Finite termination leads directly to a number of useful 
auxiliary guarantees on both the nature and quality of the 
output tessellation. Adopting the conventional terminol¬ 
ogy, surface meshes generated using the Frontal-Delaunay 
algorithm can be considered to be ‘provably-good’: 

Corollary 1. Given a smooth surface S, a radius-edge 
threshold p > 1 and positive mesh-size and surface-error 
functions h(x) > 0 , e(x) > 0 , all surface facets / G T|e i'^ 
the restricted surface tessellation T|e = Del|E(X) gener¬ 
ated by the Frontal-Delaunay algorithm satisfy: (i) p{f) < 
p, (a) h{f) < ah{xf), and (Hi) e(/) < 6 (x/). 

Proof. The Frontal-Delaunay maintains a queue of ‘bad’ 
elements, containing any surface facets / G Del |e(-^) that 
are in violation of one or more local constraints. Specif¬ 
ically, a given 2-face / G Del|E(X) is enqueued if: (i) 
p{f) > P, (ii) Kf) > ah{xf), or (iii) e(/) > e(x/). Given 
that termination is assured, it is clear that the refinement 
queue must become empty after a finite number of steps 
and that all facets in the resulting triangulation Del |e(-^) 
satisfy the requisite radius-edge ratio, element size and 
surface error constraints as a consequence. □ 

In addition to proofs of termination and convergence, 
it is also necessary to ensure the point-set X G H remains 
well-distributed throughout the refinement process. Such 
behaviour requires that no spurious short edges be intro¬ 
duced within the initial refinement phase as the algorithm 
seeks to satisfy mesh-size and surface-error constraints. 

Proposition 4. Let Del|E(X/c) be the current restricted 
surface triangulation. Given a large refinable surface facet 
f G Del|s(X/c), the minimum edge length ||eo||/c is de¬ 
creased by an 0(1) factor of in the worst-case fol¬ 
lowing the insertion of a new Type I or Type II Steiner 
vertex. 

Proof. Considering the insertion of new Type I and Type II 
vertices separately: 

(i) Recalling that the surface facet / is locally Delaunay, 
the insertion of a new Type I vertex at the centre 
of the associated surface Delaunay ball B(c/e,r/c) = 


SDB(/) is guaranteed to create edges no shorter than 
the radius of the ball of /: 


||eo||fe+i > rfc. (10) 


Here ||eo||/c+i is the length of the shortest edge cre¬ 
ated by the insertion of the new vertex c/^. Dividing 
by the current minimum edge length, and using the 
definition of the radius-edge ratio: 


l|eolU+i > rk 

lleolU “ lleolU 


P{f)- 


( 11 ) 


Clearly, the minimum edge length is decreased when 
any element with p{f) < 1 is refined. Noting that 
p{f) achieves a minimum of yVs when / is equi¬ 
lateral, it can be seen that insertion of a Type I 
Steiner vertex leads to a reduction in the minimum 
edge length by a factor of i/x/s in the worst case. 

(ii) A Type II vertex Ck positioned on a facet V C Vor(X) 
associated with the short edge eo, is equidistant to 
the vertices x^, Xj G cq. Since V is a local facet of the 
Voronoi complex, no other vertex x^ G X lies closer 
to the point C/^. 

Recalling that the Type II refinement strategy re¬ 
quires the diametric ball of the edge cq remain empty, 
a worst case reduction in ||eo|| occurs when Ck is po¬ 
sitioned at the intersection of V and the diametric 
ball of eo, forming a right triangle. Noting that V in¬ 
tersects eo at its midpoint, the minimum edge length 
is decreased by a factor of ^/V 2 in such cases. 


Overall, Type I point-placement is a limiting case, re¬ 
ducing minimum edge length by an 0 ( 1 ) factor of in 
the worst case. □ 


3.7. Quality of Approximation 

While the preceding results address the behaviour and 
performance of the underlying Delaunay-refinement algo¬ 
rithm, the approximation quality of the resulting mesh is 
also a key consideration. Using the framework introduced 
by Boissonnat and Oudot mms], a variety of geometrical 
and topological guarantees are sought. Firstly, a number 
of results from ns US] are summarised: 

Definition 5. Given a bounded domain Q enclosed by a 
surface S, the medial-axis of U is the topological closure 
of the set of centres of empty balls that touch the surface 
S at more than one point. 

Noting that the medial-axis lies at the ‘centre’ of any 
geometrical features in the domain induced by S, it is 
clear that the distance to the medial axis, denoted dM(x) 
throughout, is a measure of the local ‘thickness’ of the do¬ 
main. Using such considerations, Boissonnat and Oudot 
introduce the notion of loose e-samples to describe point- 
wise samplings X G S that induce restricted surface trian¬ 
gulations Del |e(X) exhibiting favourable geometrical and 
topological guarantees. 
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Definition 6. Let X G S be a point-wise sampling of a 
smooth surface S. Let Sy{X) C Vor(X) be the set of 
edges in the associated Voronoi complex. The sampling X 
is called a is called a loose e-sample of H if: (i) for all inter¬ 
sections X G the associated e-balls, B(c, r), 

are non-empty, such that X D B(x, edM{^)) ^ 0 for some 
e G and (ii) the restricted triangulation Del |e(-^) in¬ 
cludes at least one vertex on all connected components of 
the surface 

Boissonnat and Oudot additionally establish the fol¬ 
lowing geometrical and topological guarantees for restricted 
Delaunay tessellations Del |e(-^) induced by such loose e- 
samples: 

Proposition 5. Let X G be a loose e-sampling of a 
smooth surfaee S, with e < 0.08. The assoeiated restrieted 
surfaee triangulation Del |e(X) exhibits the following prop¬ 
erties: (i) the triangulation Del |e(-^) is homeomorphie to 
the surfaee S, (ii) for any 2-faee f G Del |e(-^) the angle 
between the faeet normal Uf and the surfaee normal 
sampled at the vertiees x e f, is 0{e), and (Hi) the one¬ 
sided Hausdorff distanee H(Del |e(X), S), sampled at the 
eentre of the surfaee Delaunay balls B(c, = SDB(/) 

is O(e^) • 

The reader is referred to [niiin] or m for additional 
details and proofs. Boissonnat and Oudot go on to show 
that restricted Delaunay refinement, consistent with the 
algorithm presented in Sectionj^ is a ‘provably-good’ tech¬ 
nique, being guaranteed to produce a loose e-sample of a 
surface H in a finite number of steps, and to generate a 
high-quality surface tessellation Del |e(X) as a result. 

The Frontal-Delaunay algorithm presented in this study 
derives from the same Delaunay-refinement framework de¬ 
veloped in [m US], and inherits much of its theoretical 
robustness as a result. Using the properties of loose e- 
samples, the surface triangulations generated using the 
Frontal-Delaunay algorithm are guaranteed to be a good 
geometrical and topological approximation to the underly¬ 
ing surface - provided the mesh size function is sufficiently 
small: 

Proposition 6. Given a smooth surfaee S and an e-eon- 
forming size funetion {^/V3)h{x.) < edM{^), the point-wise 
sampling X G S generated by the Frontal-Delaunay algo¬ 
rithm is a loose e-sample. When e < 0.08 the assoeiated 
triangulation Del|s(X) is a tight topologieal and geomet- 
rieal approximation of the surfaee S. 

Proof The Frontal-Delaunay algorithm refines any large 
surface facets, ensuring that h{f) < ah(Kf) for all / G 
Del |e(X). Rearranging the previous expression, and sub¬ 
stituting for the given mesh size constraints leads to h{f) = 
y^r < o(v^/a)e(iM(x), which simplifies to r < 0.08 dM(x) 
when e < 0.08. Recalling Definition it is clear that X is 
a loose e-sample of S, and, via Proposition]^ Del |e(-^) is 
a tight geometrical and topological approximation of S as 
a result. □ 


4. Mesh Size Functions 

The construction of high-quality mesh-size functions 
is an important aspect of both the restricted Delaunay- 
refinement and Frontal-Delaunay algorithms presented in 
Sections]^ andA good mesh size function h(x) incorpo¬ 
rates sizing constraints imposed by both the user and the 
geometry of the domain to be meshed. These contributions 
can be considered via two separate size functions, where 
hu (x) represents user-defined sizing information and hg (x) 
encapsulates geometry driven constraints. In this study, it 
is required that both hn(x) and h^(x) be piecewise linear 
functions defined on a supporting tetrahedral complex S. 
Construction of appropriate user-defined functions hn(x) 
is highly problem dependent and detailed formulations are 
not considered here. 

As per the discussions presented in Section |3.5[ it is 
known that the geometric mesh size function h^(x) must 
be sufficiently small relative to the loeal-feature-size of the 
domain, denoted Ifs(x), to ensure that the resulting sur¬ 
face triangulation Del |e(X) is both topologically and ge¬ 
ometrically correct. Specifically, to guarantee that the re¬ 
sulting tessellation is a loose e-sample, it is required that 
hg(x.) < 0.081fs(x), where Ifs(x) is the distance to the 
medial-axis of the bounded domain O. In practice, such 
constraints are typically found to be overly restrictive [19] , 
allowing for a relaxation of the coefficient e. In this study, 
e = 1/2 is used throughout. The local-feature-size can be 
expressed as the distance from any point x^ G H to the 
closest point on the medial-axis of the domain. Given that 
S is represented as an underlying triangulation V in this 
work, a discrete approximation to Ifs(x) can be calculated 
directly via the method of Amenta and Bern [2ll[3Il[38|. 
For the sake of brevity, we do not describe these methods 
in full here, but simply note that they provide an efficient 
means to estimate Ifs(x) at the vertices of the supporting 
complex S. 

Given an estimate of Ifs(x) on the boundary of 5, it 
is possible to exact a degree of user-defined control on the 
resulting size function via a Lipschitz smoothing process. 
Following the work of Persson [39|, a gradient-limited func¬ 
tion h(x) can be constructed by limiting the variation in 
h(x) over the elements of S. In this study, a scalar smooth¬ 
ing parameter g G is used to globally, and isotropically 
limit variation, such that 

h{xi) < h{xj) + g \\xi - XjW ( 12 ) 

for all vertex pairs x^,Xj G S. The gradient-limited size 
function h(x) becomes more uniform as ^ ^ 0 . 

5. Results Discussions 

The performance of the Delaunay-refinement and Frontal- 
Delaunay surface meshing algorithms presented in Sec¬ 
tions and was investigated experimentally, with both 
techniques used to mesh a series of benchmark problems. 
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Figure 4: Triangulations for the elephant, hip and bunny problems, showing output for the JGSW-FD algorithm (upper) and the CGAL-dr 
algorithm (lower). Meshes were built using: (i) uniform mesh size functions h(x) = a, with a G M+, (ii) tight constraints on element 
shape-quality, such that ^min > 30°, and (iii) uniform surface discretisation thresholds e = ( 1 / 4 ) h(x). Element counts ITIeI, and algorithm 
run-times (t) are included for each case. Normalised histograms of element area-length ratio «(/), plane-angle 0(f) and relative-length hr are 
illustrated. 
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Both the Frontal-Delaunay and Delaunay-refinement algo¬ 
rithms were implemented, allowing the performance and 
output of the two techniques to be compared side-by-side. 
Due to similarities in the overall algorithmic structure, a 
common code-base was used, with the algorithms differing 
only in the type of Steiner vertices inserted, as per the dis¬ 
cussions outlined in Section O Care was taken to ensure 
that both methods were implemented in a consistent fash¬ 
ion, allowing unbiased comparisons to be made between 
the algorithms without needing to account for systemic 
differences arising from particular implementation and/or 
design choices. Both algorithms were implemented in C++ 
and compiled as 64-bit executables. Both the Frontal- 
Delaunay and Delaunay-refinement algorithms have been 
implemented as part of the JIGSAW meshing package, cur¬ 
rently available by request from the first author. These 
implementations are referred to as JGSW-FD and JGSW- 
DR throughout the following discussions, with the suffixes 
-FD and -DR denoting ‘Frontal-Delaunay’ and ‘Delaunay- 
refinement’, respectively. 

In order to provide additional performance informa¬ 
tion, the well-known GGALMESH implementation m was 
also included in a subset of the meshing comparisons. The 
GGALMESH algorithm was sourced from version 4.6 of the 
GGAL package gni E] and was compiled as a 64-bit li¬ 
brary. The GGALMESH algorithm is referred to as GGAL- 
DR throughout the following discussions, with the sufhx 
-DR denoting ‘Delaunay-refinement’. All tests were run 
using a single core of an Intel i7 processor. Visualisation 
and post-processing was completed using matlab. 

5.1. Preliminaries 

The various Delaunay-refinement and Frontal-Delaunay 
algorithms were used to mesh a series of eight benchmark 
problems, including (i) a set of three test-cases presented 
in Figure in which the performance of the new JGSW-FD 
algorithm and the existing GGAL-dr implementation was 
compared, (ii) three additional test-cases shown in Fig¬ 
ure designed to contrast the performance of the JGSW- 
FD and JGSW-DR algorithms, and (iii) two detailed com¬ 
parative studies presented in Figures and [^designed to 
assess the impact of non-uniform mesh-sizing constraints 
on the performance of the JGSW-FD and JGSW-DR algo¬ 
rithms. 

In all test cases, a constant radius-edge ratio threshold, 
p = 1 was specified, corresponding to ^min > 30°. Ad¬ 
ditionally, non-uniform surface discretisation constraints 
were enforced, setting e(x) = /3h(x), with [3 = 1 / 4 . For 
all test problems, detailed statistics on element quality 
are presented, including histograms of element area-length 
a{f), plane-angle 0{f) and relative-length hr metrics. His¬ 
tograms further highlight the minimum and mean area- 
length metrics, the worst-case plane angle bounds, Omin 
and 6>max and the mean relative edge-length. 

Overall, both Delaunay-refinement algorithms (ggal- 
DR, jgsw-dr) and the new Frontal-Delaunay technique 
(jgsw-fd) were found to successfully mesh the full set 


of benchmark problems, satisfying all constraints. These 
results demonstrate that the new Frontal-Delaunay tech¬ 
nique is capable of generating high-quality surface trian¬ 
gulations, satisfying constraints on element shape-quality 
p(/), element size h(x/) and surface discretisation error 
e(xj), consistent with conventional Delaunay-refinement 
approaches. 


5.2. Mesh Quality Metries 

Before moving on to a detailed analysis of results, a 
number of mesh quality metrics are first introduced. 

Definition 7. Given a 2-simplex /, the plane-angle be¬ 
tween any two adjacent edges is given by: 

° WlM 

for all adjacent pairs ef,ej G /. 

Clearly, high-quality surface triangulations include a 
majority of element plane-angles clustered about 60°. 

Definition 8. Given a 2-simplex /, its area-length ratio, 
«(/)> is given by: 


«(/) = 


4\/3 Af 
3 ||e™,||2’ 


(14) 


where Af is the signed-area of / and ||erms|| is the root- 
mean-square edge length. 


The area-length ratio is a robust, scalar measure of ele¬ 
ment shape-quality, and is typically normalised to achieve 
a score of +1 for ideal elements. The area-length ratio de¬ 
creases with increasing distortion, achieving a score of -f-O 
for degenerate elements and —1 for fully inverted elements. 


Definition 9. Given a 2-simplex /, the relative-length of 
its edges is given by: 


hr{ei) 


I|e4 

h(Xe) 


(15) 


where ||ej|| is the length of the i-th edge and /i(xe) is the 
mesh size function sampled at the edge midpoint. 


The relative-length distribution hr is a measure of size- 
function conformance, expressing the ratio of actual-to- 
desired edge length for all edges e G Del A value of 

hr = 1 indicates perfect conformance. 


Definition 10. Given a distribution of element-wise plane 
angles 6>(/i), the Mean Absolute Deviation MAD (Of) is 
given by: 


MAD{ef) = l^mi)-0j\ (16) 

where 0{fi) denotes the mean value of the underlying dis¬ 
tribution 0{fi). 
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Figure 5: Triangulations for the ifp2, femur and ROCKER test problems, showing output for the JGSW-FD (upper) and JGSW-DR (lower) 
algorithms. Meshes were built using: (i) uniform mesh size functions h(x) = a, with a G M+, (ii) tight constraints on element shape-quality, 
such that ^min > 30°, and (iii) non-uniform surface discretisation functions e = ( 1 / 4 )/i(x). Element counts ITIeI, and algorithm run-times 
(t) are included for each case. Normalised histograms of element area-length ratio «(/), plane-angle 6 (f) and relative-length hr are also 
illustrated. 


(jgsw-fd): |T|e| = 87,083, t = 6.24s 


(jgsw-fd): |T|e| = 7,854, t = 0.48s 


(jgsw-fd): |T|s| = 13,023, t = 0.72s 
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MAD(^j) is a measure of the spread of the distribution 
of element angles, with smaller values indicative of distri¬ 
butions that are more tightly clustered about the mean 
value. In this study, MAD (Of) is used as a measure of av¬ 
erage element quality, with smaller values associated with 
higher quality meshes. 

5.3. A Comparison of JGSW-FD and CGAL-dr 

The performance of the JGSW-FD and GGAL-dr algo¬ 
rithms was assessed using a set of three test-cases. The 
ELEPHANT, HIP and BUNNY geometries were meshed us¬ 
ing uniform mesh-size constraints, with a small constant 
value, h(x) = a, imposed globally, where a was chosen 
to be approximately 2% of the mean bounding-box di¬ 
mension associated with each model. Due to differences 
in the way that mesh-size constraints are interpreted by 
the two meshing packages, the mesh-size value selected 
for the GGAL-DR algorithm was reduced by a factor of 
4 / 3 . This reduction accounts for the fact that in GGAL-DR, 
mesh-size constraints are imposed with respect to the ele¬ 
ment circumradii, whereas the JGSW-FD algorithm treats 
the mesh-size function in terms of element edge length. 
Both algorithms were found to produce meshes incorporat¬ 
ing consistent mean edge-lengths based on these modified 
mesh-size values. 

Overall, the results shown in Figure demonstrate 
that both the JGSW-FD and GGAL-DR algorithms gener¬ 
ate high-quality surface meshes for all test cases - satis¬ 
fying the required element shape-quality, mesh-size and 
surface discretisation error thresholds. Focusing on the 
distribution of element shape-quality explicitly, it is clear 
that the JGSW-FD algorithm achieves significantly better 
results, producing meshes with higher mean area-length 
ratios in all cases. In fact, with mean area-length ratios 
of d(/) 0.98 (compared to d(/) 0.92 for GGAL-dr), 

it is clear that the vast majority of elements generated by 
the JGSW-FD algorithm are of near-ideal shape. This im¬ 
provement in mean mesh-quality is also evident in the dis¬ 
tribution of element plane-angles, with the JGSW-FD ap¬ 
proach generating histograms clustered more tightly about 
60°. These results are confirmed by an analysis of the 
mean-absolute-deviation in angle, MAD (Of) - included in 
Table showing that the spread of the plane-angle dis¬ 
tribution is reduced by a factor of approximately three 
through use of the JGSW-FD algorithm. Finally, analy¬ 
sis of the relative-length distributions shows that meshes 
generated using the JGSW-FD algorithm tightly conform 
to the imposed mesh size constraints, with a tight cluster¬ 
ing of hr about 1. In contrast, output generated using the 
GGAL-DR scheme is seen to incorporate significant sizing 
‘error’, typified by broad distributions of relative-length, 
straddling A — 1. 

In terms of computational expense, meshes generated 
by the JGSW-FD and GGAL-DR algorithms are shown to 
be very similar in size. Additionally, despite the additional 
work required to support the off-centre refinement scheme. 


the JGSW-FD algorithm appears to be slightly more effi¬ 
cient than GGAL-DR, though the difference in run-time 
performance is not significant. 

5 . 4 . A Comparison of JGSW-FD and JGSW-DR 

Consistent with the analysis presented in Section |5.3[ 
the performance of the JGSW-FD and JGSW-DR algorithms 
was next assessed. While GGAL-DR and JGSW-DR are rep¬ 
resentative of the same class of meshing algorithms^ the 
benchmarks included in this section allow for an unbi¬ 
ased comparison of the Frontal-Delaunay and Delaunay- 
refinement approaches, with both the JGSW-FD and JGSW- 
DR algorithms benefiting from the same set of implemen¬ 
tation design and optimisation decisions. 

Consistent with the previous analysis, the performance 
of the JGSW-FD and JGSW-DR algorithms was assessed 
using a set of three test-cases. The ifp2, femur and 
ROGKER geometries shown in Figure were meshed us¬ 
ing uniform mesh-size constraints, with a small constant 
value, h(x) = a, imposed globally, where a was chosen to 
be approximately 2% of the mean bounding-box dimen¬ 
sion associated with each model. An analysis of the dis¬ 
tributions of element area-length, plane-angle and relative 
edge-length measures confirms much the behaviour noted 
in Section |5.3| - in terms of overall mesh-quality and size- 
conformance, the JGSW-FD algorithm is a clear winner. 

An analysis of computational expense leads to different 
conclusions, with the JGSW-DR algorithm seen to produce 
meshes that are slightly and yet consistently larger than 
those generated using JGSW-FD. The JGSW-DR scheme 
also appears to be somewhat faster than the JGSW-FD al¬ 
gorithm, typically requiring only approximately 70% of the 
run-time imposed by JGSW-FD. Such behaviour is not un¬ 
expected, with the additional work associated with the off- 
centre point-placement scheme employed in JGSW-FD ex¬ 
pected to result in lower overall performance. The search 
for additional efficiency in the JGSW-FD algorithm is an 
interesting topic for future research. 

5.5. Non-uniform Mesh-Size Constraints 

A sequence of meshes for the BIMBA and KISS prob¬ 
lems are presented in Figures and examining the im¬ 
pact of non-uniform mesh-size constraints on algorithm 
performance. A set of meshes were generated for both 
test cases using low, medium and high resolution settings, 
where the associated mesh size functions h(x) were built 
using increasingly stringent gradient limits, such that gi G 
{3/10, 2/10,1/10} respectively. Analysis of Figures and 
show that both the JGSW-FD and JGSW-DR algorithms 
generate high-quality meshes for all test cases - satisfying 
the required element quality, size and surface error thresh¬ 
olds. Analysis of the area-length and plane-angle distribu¬ 
tions again shows that the JGSW-FD algorithm consistently 
outperforms the JGSW-DR scheme, generating meshes with 
higher mean area-length ratios in all cases. Furthermore, 
it is evident that the quality of meshes generated using the 
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Figure 6: Size-driven meshing for the BIMBA problem, showing output for the JGSW-FD (upper) and JGSW-DR (lower) algorithms. Meshes 
were built using: (i) graded, ^f-Lipschitz smooth mesh size functions h(x), with gi G {3/io, 2 / 10 , ^/lo} from left to right, (ii) tight constraints 
on element shape-quality, such that ^min ^ 30°, and (iii) non-uniform surface discretisation functions e = ( 1 / 4 ) h(x). Element counts ITIeI? 
and algorithm run-times (t) are included for each case. Normalised histograms of element area-length ratio «(/), plane-angle 0{f) and 
relative-length hr are illustrated. 


(jgsw-fd): |T|s| = 41,404, t = 3.10s 


(jgsw-fd): |T|e| = 48,860, t = 3.70s 


(jgsw-fd): |T|e| = 69,422, t = 5.35s 
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Figure 7: Size-driven meshing for the KISS problem, showing output for the JGSW-FD (upper) and JGSW-DR (lower) algorithms. Meshes were 
built using: (i) graded, ^f-Lipschitz smooth mesh size functions h(x), with gi G {3/io, 2 / 10 , 1 / 10 } from left to right, (ii) tight constraints on 
element shape-quality, such that ^min > 30°, and (iii) non-uniform surface discretisation functions e = ( 1 / 4 ) h(x). Element counts ITIeI, 
and algorithm run-times (t) are included for each case. Normalised histograms of element area-length ratio «(/), plane-angle 0{f) and 
relative-length hr are illustrated. 
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Table 1: A comparison of the mean absolute deviation in the element angle distributions for the various combinations of refinement 

algorithms and benchmark problems considered in this study. Smaller values correspond to improved mean mesh quality. CGAL-dr denotes the 
Delaunay-refinement algorithm included in the CGAL package, and JGSW-dr/jgsw-fd denote the Delaunay-refinement and Frontal-Delaunay 
algorithms implemented in the current study. The rightmost column indicates the reduction in absolute deviation achieved using the Frontal- 
Delaunay algorithm. 
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GGAL-DR 

JGSW-DR 

JGSW-FD 

FAGTOR 

ELEPHANT 

- 

11.96° 

- 

4.35° 

2.8 

HIP 

- 

11.90° 

- 

3.89° 

3.1 

BUNNY 

- 

12.02° 

- 

4.10° 

2.9 

ifp2 

- 

- 

10.72° 

5.14° 

2.1 

FEMUR 

- 

- 

10.93° 

3.85° 

2.8 

ROGKER 

- 

- 

10.69° 

3.92° 

2.7 

BIMBA 

3/l0 

- 

11.25° 

7.83° 

1.4 


2/l0 

- 

10.99° 

7.04° 

1.6 


i/io 

- 

10.74° 

5.68° 

1.9 

KISS 

3/l0 

- 

11.37° 

8.15° 

1.4 


2/l0 

- 

10.98° 

7.16° 

1.5 


i/io 

- 

10.74° 

5.74° 

1.9 


Frontal-Delaunay algorithm improves with increased uni¬ 
formity, via ^ 0 , as indicated by the increase in mean 
a(/) and the narrowing of the 0{f) distribution about 
60.0°. Results for the mean absolute deviation in the ele¬ 
ment angle distribution, MAD (Of) - included in Table 
shows that use of the JGSW-FD algorithm results in a re¬ 
duction in the spread of the element angle distributions by 
a factor ranging from approximately 1.5 to 2 . 0 . 

In contrast, similar analysis shows that output gener¬ 
ated via the JGSW-DR scheme to be essentially indepen¬ 
dent of sizing uniformity, with distributions of a{f) and 
0{f) showing little variation with h(x). Visually, the en¬ 
hanced quality of the meshes generated using the Frontal- 
Delaunay algorithm is again evident, with all test cases 
showing a marked increase in both smoothness and sub¬ 
structure. Given the tight bounds on radius-edge ratios in 
these tests, the Delaunay-refinement algorithm appears to 
achieve a maximum mean area-length ratio a{f) 0.92. 
The Frontal-Delaunay algorithm is seen to generate consis¬ 
tently better output, with mean area-length metrics reach¬ 
ing d(/) ^ 0.97 at high resolution settings. 

Consistent with previous results, analysis of the dis¬ 
tribution of relative-length ratios show that the JGSW-FD 
algorithm tightly conforms to the imposed non-uniform 
sizing constraints, with hr tightly clustered about 1. Fur¬ 
thermore, like the element shape-quality results discussed 
previously, it can be seen that this distribution improves 
as ^ 0 , as indicated by the narrowing of the distribu¬ 
tion about = 1. In contrast, the Delaunay-refinement 
scheme is again shown to be consistently associated with 
significant sizing error, as illustrated by broad distribu¬ 
tions of relative-length straddling — I. Such behaviour 
appears to be largely independent of the characteristics of 
the imposed mesh size constraints. 


6. Conclusions 

A new Frontal-Delaunay algorithm has been developed 
to triangulate smooth surfaces embedded in The new 
algorithm is based on the so-called restricted Delaunay 
paradigm, in which a surface mesh Del conforming 

to an underlying surface description S, is constructed as a 
subset of a full-dimensional Delaunay tessellation Del(X). 
The new restricted Frontal-Delaunay algorithm is based 
on a generalisation of the work of Rebay [28] and Ungor 
[3^ for planar problems, in which generalised ‘off-centre’ 
Steiner vertices are inserted along edges in the Voronoi 
diagram. This work has been extended to support sur¬ 
face meshing operations through the development of a new 
point-placement strategy that positions vertices on facets 
in the underlying Voronoi complex. This new scheme 
allows for the insertion of both size- and shape-optimal 
Steiner vertices, leading to a hybrid approach that com¬ 
bines many of the advantages of conventional advancing- 
front and Delaunay-refinement techniques. A series of 
comparative experimental studies confirm the effectiveness 
of this new approach in practice, demonstrating that an 
improvement in mesh-quality is typically achieved when 
compared to conventional Delaunay-refinement schemes. 
Importantly, it has also been demonstrated that the new 
Frontal-Delaunay algorithm satisfies the same set of con¬ 
straints as conventional restricted Delaunay-refinement ap¬ 
proaches, adhering to limits on element radius-edge ratios, 
edge length and surface discretisation error. Results show 
that the new algorithm is an effective hybridisation of ex¬ 
isting mesh generation techniques, combining the high el¬ 
ement quality and mesh sub-structure of advancing-front 
techniques with the theoretical guarantees of Delaunay- 
refinement schemes. It is expected that applications that 
place a premium on mesh-quality, including problems in 
computational fluid dynamics and/or structural analysis, 
may benefit from the new Frontal-Delaunay technique. Fu- 
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ture work should focus on support for an extended class 
of surface definitions, including domains containing sharp 
features. 

Acknowledgements 

This work was carried out at the University of Sydney 
with the support of an Australian Postgraduate Award. 
The authors also wish to thank the anonymous reviewers 
for their helpful comments and feedback. 

References 

References 

[1] D. Engwirda, D. Ivers, Face-centred Voronoi refinement for sur¬ 
face mesh generation, Procedia Engineering 82 (2014) 8-20. 

[2] M. Bern, D. Eppstein, J. Gilbert, Provably Good Mesh Gener¬ 
ation, J. Gomput. Syst. Sci 48 (1990) 231-241. 

[3] S. A. Mitchell, S. A. Vavasis, Quality Mesh Generation in 
Three-dimensions, in: Symposium on Gomputational Geome¬ 
try, 1992, pp. 212-221. 

URL http: / / dblp. uni-trier. de/db/conf / compgeom/ 

compgeom92.html#MitchellV92 

[4] J. Peraire, J. Peiro, K. Morgan, Advancing-Front Grid Gener¬ 
ation in: J. F. Thompson, B. K. Soni, N. P. Weatherill (Eds.), 
Handbook of Grid Generation, Taylor & Francis, 1998. 

URL http: //books. google. com. au/books?id=ImaDT6i jKq4C 

[5] J. Schoberl, NETGEN: An Advancing Front 2D/3D Mesh Gen¬ 
erator based on Abstract Rules Gomputing and Visualization 
in Science 1 (1) (1997) 41-52. doi:10.1007/s007910050004 
URL http: //dx. doi. org/10.1007/s007910050004 

[6] D. Rypl, Approaches to Discretization of 3D Surfaces, in: GTU 
Reports, Vol. 7, GTU Publishing House, Prague, Gzech Repub¬ 
lic, 2003. 

[7] J. Schreiner, G. E. Scheidegger, S. Fleishman, G. T. Silva, Direct 
(Re)Meshing for Efficient Surface Processing Gomputer Graph¬ 
ics Forum 25 (3) (2006) 527-536. doi: 10.1111/j . 1467-8659. 
2006.00972.x 

URL http: //dx. doi. org/10.1111/j . 1467-8659.2006.00972. x 

[8] L. P. Ghew, Guaranteed-quality Triangular Meshes, Tech, rep., 
Gornell University, Department of Gomputer Science, Ithaca, 
New York (1989). 

[9] J. Ruppert, A New and Simple Algorithm for Quality 2- 
dimensional Mesh Generation in: Proceedings of the fourth 
annual AGM-SIAM Symposium on Discrete algorithms, SODA 
’93, Society for Industrial and Applied Mathematics, Philadel¬ 
phia, PA, USA, 1993, pp. 83-92. 

URL http: //dl. acm. org/citation. cfm?id=313559.313615 

[10] J. Ruppert, A Delaunay Refinement Algorithm for Quality 2- 
Dimensional Mesh Generation Journal of Algorithms 18 (3) 
(1995) 548 - 585. doi:http://dx.doi.org/10.1006/jagm.1995. 
1021 

URL http://www.sciencedirect.com/science/article/pii/ 
S0196677485710218 

[11] J. R. Shewchuk, Delaunay Refinement Mesh Generation, Ph.D. 
thesis, Pittsburg, Pennsylvania (1997). 

[12] J. R. Shewchuk, Tetrahedral Mesh Generation by Delaunay Re¬ 
finement in: Proceedings of the fourteenth annual symposium 
on Gomputational geometry, SGG ’98, AGM, New York, NY, 
USA, 1998, pp. 86-95. doi:10.1145/276884.276894 

URL http: //doi. acm. org/10.1145/276884.276894 

[13] S. Gheng, T. Dey, Quality Meshing with Weighted Delaunay 

Refinement SIAM Journal on Gomputing 33 (1) (2003) 
69-93. arXiv:http://epubs.siam.org/doi/pdf/lO.1137/ 

S0097539703418808 doi:10.1137/S0097539703418808 

URL http://epubs.siam.org/doi/abs/10.1137/ 

S0097539703418808 


[14] S. W. Gheng, T. K. Dey, E. A. Ramos, Delaunay Refine¬ 
ment for Piecewise Smooth Gomplexes Discrete & Gompu¬ 
tational Geometry 43 (1) (2010) 121-166. doi: 10.1007/ 
S00454-008-9109-3 

URL http: //dx. doi. org/10.1007/s00454-008-9109-3 

[15] J. D. Boissonnat, S. Oudot, Provably Good Surface Sampling 
and Approximation, in: AGM International Gonference Pro¬ 
ceeding Series, Vol. 43, 2003, pp. 9-18. 

[16] J. D. Boissonnat, S. Oudot, Provably Good Sampling and Mesh¬ 
ing of Surfaces, Graphical Models 67 (5) (2005) 405-451. 

[17] G. Jamin, P. Alliez, M. Yvinec, J. D. Boissonnat, GGALmesh: 
A Generic Framework for Delaunay Mesh Generation, Tech, 
rep., INRIA (2013). 

[18] B. Delaunay, Sur la sphere vide, Izv. Akad. Nauk SSSR, Otde- 
lenie Matematicheskii i Estestvennyka Nauk 7 (793-800) (1934) 
1-2. 

[19] S. W. Gheng, T. K. Dey, J. R. Shewchuk, Delauay Mesh Gen¬ 
eration, Taylor & Francis, New York, 2013. 

[20] H. Edelsbrunner, N. R. Shah, Triangulating Topological Spaces, 
International Journal of Gomputational Geometry & Applica¬ 
tions 7 (04) (1997) 365-378. 

[21] N. Amenta, M. Bern, Surface Reconstruction by Voronoi Filter¬ 
ing, Discrete & Gomputational Geometry 22 (4) (1999) 481-504. 

[22] S. W. Gheng, T. K. Dey, J. A. Levine, A Practical Delau¬ 
nay Meshing Algorithm for a Large Glass of Domains in: 
M. Brewer, D. Marcum (Eds.), Proceedings of the 16th Interna¬ 
tional Meshing Roundtable, Springer Berlin Heidelberg, 2008, 
pp. 477-494. doi:10.1007/978-3-540-75103-8\_27 

URL http: //dx. doi. org/10.1007/978-3-540-75103-8_27 

[23] S. W. Gheng, T. K. Dey, E. A. Ramos, T. Ray, Sampling and 
Meshing a Surface with Guaranteed Topology and Geometry, 
SIAM journal on computing 37 (4) (2007) 1199-1227. 

[24] T. K. Dey, J. A. Levine, Delaunay meshing of piecewise 
smooth complexes without expensive predicates. Algorithms 
2 (4) (2009) 1327-1349. 

[25] P. A. Foteinos, A. N. Ghernikov, N. P. Ghrisochoides, Guaran¬ 
teed quality tetrahedral Delaunay meshing for medical images, 
Gomputational Geometry: Theory and Applications 47 (4) 
(2014) 539-562. 

[26] A. Bowyer, Gomputing Dirichlet Tessellations, The Gomputer 
Journal 24 (2) (1981) 162-166. 

[27] H. Erten, A. Ungor, Quality Triangulations with Locally Op¬ 
timal Steiner Points, SIAM J. Sci. Gomp. 31 (3) (2009) 2103- 
2130. 

[28] S. Rebay, Efficient Unstructured Mesh Generation by Means of 
Delaunay Triangulation and Bowyer-Watson Algorithm Jour¬ 
nal of Gomputational Physics 106 (1) (1993) 125 - 138. doi: 
http://dx.doi.org/10.1006/j cph.1993.1097 

URL http: //WWW. sciencedirect. com/science/article/pii/ 
S0021999183710971 

[29] D. J. Mavriplis, An Advancing Front Delaunay Triangulation 
Algorithm Designed for Robustness Journal of Gomputational 
Physics 117 (1) (1995) 90 - 101. doi:http://dx.doi.org/10. 
1006/jcph.1995.1047 

URL http: //www .sciencedirect.com/science/article/pii/ 
S0021999185710479 

[30] P. J. Frey, H. Borouchaki, P. L. George, 3D Delaunay Mesh Gen¬ 
eration Goupled with an Advancing-Front Approach Gomputer 
Methods in Applied Mechanics and Engineering 157 (12) (1998) 
115 - 131. doi:http://dx.doi.org/10.1016/30045-7825(97) 
00222-3 

URL http: //www .sciencedirect.com/science/article/pii/ 
S0045782597002223 

[31] J. F. Remade, F. Henrotte, T. Garrier-Baudouin, E. Bechet, 
E. Marchandise, G. Geuzaine, T. Mouton, A Frontal-Delaunay 
Quad-Mesh Generator using the L°°-norm International Jour¬ 
nal for Numerical Methods in Engineering 94 (5) (2013) 494- 
512. doi:10.1002/nme.4458 

URL http: //dx. doi. org/10.1002/nme. 4458 

[32] D. Rypl, Sequential and Parallel Generation of Unstructured 
3D Meshes, Ph.D. thesis, GTU in Prague (1998). 


18 


[33] J. Schreiner, C. Scheiclegger, C. Silva, High-Quality Extraction 
of Isosurfaces from Regular and Irregular Grids, Visualization 
and Computer Graphics, IEEE Transactions on 12 (5) (2006) 
1205-1212. doi:10.1109/TVCG.2006.149 

[34] A. Ungor, Off-centers: A New Type of Steiner Points for 
Computing Size-optimal Guaranteed-quality Delaunay Trian¬ 
gulations, Computational Geometry: Theory and Applications 
42 (2) (2009) 109-118. 

[35] A. N. Chernikov, N. P. Chrisochoides, Generalized insertion re¬ 
gion guides for Delaunay mesh refinement, SIAM Journal on 
Scientific Computing 34 (3) (2012) A1333-A1350. 

[36] P. A. Foteinos, A. N. Chernikov, N. P. Chrisochoides, Fully 
generalized two-dimensional constrained Delaunay mesh refine¬ 
ment, SIAM Journal on Scientific Computing 32 (5) (2010) 
2659-2686. 

[37] N. Amenta, S. Choi, R. K. Kolluri, The power crust, unions of 
balls, and the medial axis transform. Computational Geometry 
19 (2) (2001) 127-153. 

[38] T. K. Dey, W. Zhao, Approximating the medial axis from the 
Voronoi diagram with a convergence guarantee, Algorithmica 
38 (1) (2004) 179-200. 

[39] P.-O. Persson, Mesh Size Functions for Implicit Geometries 
and PDE-based Gradient Limiting Engineering with Comput¬ 
ers 22 (2) (2006) 95-109. doi:10.1007/s00366-006-0014-1 
URL http: //dx. doi. org/10.1007/s00366-006-0014-1 

[40] C. Jamin, S. Pion, M. Teillaud, 3D Triangulations in: CGAL 
User and Reference Manual, 4.6.1 Edition, CGAL Editorial 
Board, 2015. 

URL http: //doc . cgal. org/4.6.1/Manual/packages. html# 

PkgTriangulation3Summary 

[41] C. Jamin, S. Pion, M. Teillaud, 3D Triangulation Data Struc¬ 
ture in: CGAL User and Reference Manual, 4.6.1 Edition, 
CGAL Editorial Board, 2015. 

URL http: //doc . cgal. org/4.6.1/Manual/packages. html# 

PkgTDS3Suinmary 


19 


